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1.0  INTRODUCTION 


This  report  summarizes  the  research  accomplishments  achieved  under  Grant  AFOSR 
F49620-94-1 -01 68,  “Large  Eddy  and  Unstructured  Grid  Methods  for  Heat  Transfer 
Applications  in  Propulsion  Systems.”  The  objective  of  the  research  was  to  address  some 
major  technological  needs  that  were  and  still  are  limiting  progress  in  the  development  of 
improved  design  and  predicting  capability  for  modem  high  performance  propulsion 
systems.  The  most  urgent  of  these  is  the  need  for  better  methods  for  predicting  the 
thermal  features  of  turbulent  flows  that  are  technologically  important.  Examples  of  such 
flows  occur  in  turbine  and  combustor  components  of  gas  turbine  engines.  Existing 
methods  were  considered  inadequate  for  the  accurate  prediction  of  turbulent  or 
transitional  flow  in  configurations  of  practical  interest. 

Considering  the  recent  and  planned  advances  in  computer  hardware,  it  was  believed 
that  the  time  was  right  for  increased  research  in  direct  (DNS)  and  large  eddy  simulations 
(LES),  since  more  and  more  flows  of  practical  interest  were  coming  within  reach  of 
current  and  planned  computers.  Extrapolating  computer  hardware  advances  into  the 
future,  it  appears  only  a  matter  of  time  before  LES  becomes  a  viable  alternative  for 
computing  many  turbulent  flows.  It  was  believed  to  be  important  to  bring  this  emerging 
technology  to  bear  on  urgent  problems  in  aeropropulsion  systems.  The  main  goal  of  the 
research,  then,  was  to  expand  LES  methodology  to  a  broader  class  of  flows  with  heat 
transfer  to  ultimately  contribute  to  the  solution  of  some  of  the  urgent  problems  that  are 
limiting  performance  of  aeropropulsion  systems. 
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To  obtain  flow  information  computationally,  one  has  three  possible  ways  to  proceed. 
Most  common  in  dealing  with  applied  problems  are  the  approaches  based  on  the 
Reynolds- averaged  Navier-Stokes  (RANS)  equations.  Here  turbulence  modeling  plays  a 
major  role  and  results  are  uncertain  or  unsatisfactory  for  many  complex  flows.  It  is 
widely  accepted  that  improvements  are  essential.  The  second  alternative  is  direct 
numerical  simulation  (DNS).  It  is  well  established  now  that  the  three-dimensional 
Navier-Stokes  equations  along  with  appropriate  forms  of  the  continuity  and  energy 
equations  govern  turbulent  flow,  but  limited  computer  resources  prohibit  the  resolution  of 
flows  characteristic  of  most  applications  by  DNS.  This  situation  is  destined  to  change 
eventually.  As  computer  hardware  and  algorithms  improve,  the  frontier  will  be 
continuously  pushed  back,  allowing  flows  of  increasing  practical  interest  to  be  computed 
by  direct  numerical  simulations  (DNS).  The  third  alternative  is  LES,  in  which  modeling 
is  required  only  for  the  smallest  and  most  universal  scales.  This  is  another  approach  to 
simulating  turbulence  that  will  benefit  increasingly  with  improvements  in  computer 
hardware.  In  fact,  it  is  believed  that  the  time  is  right  for  increased  research  in  this  area  as 
more  and  more  practical  flow  conditions  come  within  reach  of  current  and  planned 
computers.  The  motivation  here,  in  addition  to  responding  to  a  need  for  basic 
understanding  of  turbulent  flows,  is  the  desire  to  move  the  LES  capability  toward 
geometries  and  conditions  more  in  line  with  those  occurring  in  critical  application  areas. 

Applying  DNS  and  LES  to  more  and  more  realistic  configurations  involves  more 
than  simply  waiting  for  the  computing  power  to  become  available  and  then  applying  a 
well-established  algorithm.  Most  of  the  configurations  addressed  by  researchers  to  date 
have  allowed  the  use  of  periodic  boundary  conditions  at  inflow  and  outflow.  Such 
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conditions  free  the  researcher  from  the  task  of  specifying  specific  detailed  conditions  at 
both  inflow  and  outflow.  Effective  and  efficient  ways  of  dealing  with  more  general 
inflow  and  outflow  conditions  for  DNS  and  LES  simulations  are  generally  not  known  or 
have  likely  not  been  optimized  to  make  the  best  use  of  computational  resources. 
Modeling  the  subgrid-scale  (SGS)  transport  in  LES  is  another  issue  that  needs  further 
development  and  understanding  as  more  and  more  complex  flows  are  addressed. 
Because  most  of  the  flows  computed  by  LES  to  date  have  been  geometrically  simple  and 
at  low  Reynolds  numbers,  the  limitations  of  existing  subgrid-scale  models  have  not  been 
severely  tested. 

Several  facets  of  LES  have  been  explored  under  this  research  program.  The  results 
of  interest  to  the  technical  community  can  be  listed  in  the  following  categories; 

•  Observations  and  developments  on  algorithms  for  LES. 

•  New  results  for  incompressible  flow  in  a  square  duct. 

•  Advances  in  the  use  of  unstructured  and  zonal  embedded  grids. 

•  New  results  for  channel  flows  with  heat  transfer  and  significant  property 
variations. 

•  Some  innovations  on  modeling  inflow  and  outflow  conditions  for  LES. 

The  next  section  will  summarize  the  results  from  this  research  program  citing 
publications  where  additional  details  can  be  found.  Copies  of  dissertations  can  be 
obtained  from  UMI  Dissertation  Services,  Ann  Arbor,  Michigan.  Copies  of  computer 
codes  can  be  obtained  from  Richard  Fletcher,  Department  of  Mechanical 
Engineering,  Iowa  State  University. 
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2.0  SUMMARY  OF  RESULTS 

2.1  Overview 

Major  results  and  innovations; 

•  A  finite  volume  solution  scheme  for  the  compressible  Navier-Stokes  equations 
has  been  developed  for  large  eddy  simulation  of  turbulent  flow  using  Mach 
number  preconditioning  to  enable  efficient  solutions  of  variable  property  gas 
flows  at  low  Mach  numbers. 

•  Two  promising  procedures  have  been  found  to  eliminate  even-odd  decoupling  of 
solutions  that  often  occur  with  the  use  of  colocated  grids  and  central  difference 
numerical  schemes.  Grid  level  decoupling  (wiggles)  was  eliminated  without 
deterioration  of  the  accuracy  of  the  turbulence  statistics. 

•  A  large  eddy  simulation  of  a  channel  flow  with  one  wall  heated  and  the 
opposing  wall  cooled  at  a  temperature  ratio  of  3  has  revealed  several  interesting 
features  of  variable  property  flows  including  confirmation  of  density  rms 
percentage  fluctuations  of  9%  for  a  flow  at  M  =  0.01. 

•  For  the  first  time,  an  LES  of  channel  flow  with  uniform  heat  flux  and  significant 
property  variations  for  both  heating  and  cooling  has  been  obtained  allowing 
comparisons  with  experimentally-based  correlations  in  the  literature  and  a 
determination  of  structural  changes  induced  by  strong  property  variations. 

•  Uniform  heat  flux  studies  revealed  that  distributions  of  mean  velocities, 
temperature,  and  rms  fluctuations  of  velocities  and  temperature  were  strongly 
influenced  by  property  variations  but  could  be  nearly  collapsed  to  constant 
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property  distributions  through  the  use  of  semi-local  coordinates  as  explained 
below. 

•  Work  has  been  initiated  toward  the  use  of  unstructured  grids  in  LES  and  the 
performance  of  compressible  LES  formulations  on  regular  and  unstructured 
grids  using  both  a  dynamic  subgrid-scale  model  and  the  Smagorinsky  model 
have  been  compared  for  the  decay  of  isotropic  turbulence.  Most  recently,  the 
use  of  hexahedral  control  volumes  in  zonal  embedded  manner  has  shown  much 
promise  as  a  means  of  economically  achieving  simulations  at  high  Reynolds 
numbers.  The  scheme  allows  for  small  volumes  near  solid  boundaries  but 
permits  coarsening  in  all  three  directions.  It  is  possible  to  obtain  good  accuracy 
with  a  relatively  small  number  of  cells.  It  is  believed  that  a  major  advance  in  the 
simulation  of  high  Reynolds  number  flows  is  in  progress. 

2.2  Numerical  Strategy 

The  numerical  formulation  used  in  the  research  was  based  on  a  finite  volume 

r 

discretization  of  the  Favre-filtered  compressible  Navier-Stokes  equations.  The  approach 
is  thought  to  be  somewhat  novel  in  three  respects.  First,  there  have  been  only  a  few 
studies  reported  in  the  literature  to  date  in  which  a  fully  coupled,  compressible 
formulation  has  been  used  for  either  DNS  or  LES.  Second,  an  all-speed  strategy  is  being 
followed  in  the  numerical  simulations  which  permits  the  same  general  methodology  to  be 
applied  to  fully  incompressible  flows,  compressible  flows  at  very  low  Mach  numbers 
where  effects  of  property  variations  need  to  be  accounted  for,  and  at  transonic  and 
supersonic  speeds.  Basically,  this  flexibility  is  achieved  by  employing  a  coupled, 
compressible  formulation  with  low  Mach  number  preconditioning  [1].  With 
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preconditioning,  the  numerical  stiffness  and  slow  convergence  associated  with  traditional 
compressible  formulations  are  avoided  and,  in  fact,  the  convergence  rate  becomes 
virtually  independent  of  Mach  number  in  the  low  speed  regime.  However,  for  purely 
incompressible  simulations,  the  scheme  can  be  further  streamlined  by  specializing  it  to  a 
pseudo-compressibility  formulation.  This  avoids  the  necessity  of  solving  the  energy 
equation  for  those  flows  that  are  essentially  isothermal.  To  this  investigator’s  knowledge, 
no  other  DNS  or  LES  scheme  has  this  generality.  Third,  effects  of  heat  transfer  were 
considered.  Although  some  heat  transfer  results  have  appeared  in  the  literature,  available 
results  are  very  limited,  especially  for  separated  flows,  and  usually  do  not  take  into 
account  effects  of  variations  in  fluid  properties. 

Within  this  framework,  several  schemes  that  differ  in  detail  have  been  employed. 
The  variations  arose  due  to  searches  for  increased  accuracy,  economy,  and  the  ability  to 
accommodate  geometric  complexity.  This  has  led  to  evaluations  of  upwind  and  central 
schemes  in  the  search  for  accuracy.  Variations  suitable  for  vector  machines  such  as  the 
Cray  C90  and  parallel  systems  such  as  the  IBM  SP2  and  SGI  Origin  2000  were 
considered  in  the  search  for  economy.  Use  of  unstructured  and  zonal  embedded  grids 
was  motivated  by  interest  in  computational  economy  and  the  long-term  need  to  deal  with 
complex  geometry. 

Further  information  on  the  numerical  formulations  employed  and  lessons  learned 
will  be  given  in  the  following  section. 

2.3  Observations  on  Discretization  Methods 
The  first  lesson  learned  was  that  performance  on  traditional  laminar  flows  in  not  a 
reliable  indication  of  performance  to  be  expected  for  use  in  LES.  One  reason  for  this  is 
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that  achieving  accurate  LES  results  requires  that  very  many  length  and  time  scales  be 
accurately  resolved  whereas  most  laminar  flow  problems  are  characterized  by  relatively 
few  critical  scales.  This  observation  was  pointed  out  in  dissertation  research  completed 
by  Wang  [2,3],  where  it  was  noted  that  calculations  performed  on  both  staggered  and 
colocated  grids  gave  identical  results  for  several  laminar  test  cases  but  the  staggered  grid 
results  were  better  and  the  scheme  more  robust  when  applied  for  LES.  Further,  an 
upwind  scheme  employing  third-order  upwinding  for  advective  terms  and  fourth-order 
central  viscous  discretization  gave  very  accurate  results  in  laminar  test  cases  but  often 
failed  to  out-performed  a  second-order  central  representation  for  LES  applications.  The 
upwind  scheme  was  found  to  be  the  most  robust  scheme,  however.  Central  difference 
schemes,  especially  on  colocated  grids,  were  observed  to  fail  due  to  apparent  instability 
for  some  simulations  early  in  the  research  program. 

For  completeness,  it  should  be  mentioned  that  more  recent  studies  [4]  have 
demonstrated  that  some  of  the  stability  problems  encountered  by  Wang  was  likely  due  to 
choices  made  in  interpolating  fluxes  to  cell  faces.  Wang  [2,3]  interpolated  the  fluxes 
themselves,  but  Chidambaram  [4]  was  able  to  show  that  constmction  of  cell  face  fluxes 
from  the  interpolated  primitive  variables  was  more  robust  than  interpolating  the  fluxes. 
When  the  central  difference  schemes  of  Wang  [2,3]  were  modified  to  compute  the  fluxes 
from  interpolated  primitive  variables,  robustness  appeared  to  be  restored  and  a  failure  due 
to  stability  on  regular  grids  has  not  been  observed  to  date.  Chidambaram  [4]  discusses 
the  probable  cause  for  this  behavior. 

From  time  to  time  when  simulating  low-speed  flows  on  a  colocated  grid,  decoupling 
of  the  solution  can  be  observed  giving  rise  to  grid-frequency  oscillations  or  “wiggles”  in 
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the  solution,  particularly  the  pressure.  This  is  a  very  well  known  effect  encountered  by 
many  investigators.  The  cause  is  associated  with  the  physical  requirement  to  represent 
the  pressure  gradient  with  a  central  difference  at  low  speeds.  This  effect  can  be 
eliminated  with  the  use  of  a  staggered  grid;  however,  colocated  grids  are  much  more 
convenient  to  deal  with,  especially  when  computing  flows  in  complex  geometries  or 
when  using  unstructured  grids.  In  our  research  the  wiggles  were  not  always  a  problem, 
but  tended  to  occur  when  dealing  with  a  complex  geometry  like  the  ribbed  channel  or 
when  simulating  flows  at  high  Reynolds  numbers.  Use  of  traditional  artificial  dissipation 
to  remove  such  oscillations  typically  causes  significant  deterioration  in  the  accuracy  of 
LES  calculations. 

Two  promising  solutions  to  the  decoupling  were  developed  in  the  research. 
Chidambaram  [4]  developed  a  momentum  interpolation  correction  procedure  that 
effectively  coupled  velocities  and  pressure.  A  similar  procedure  has  been  developed  by 
Rhie  and  Chow  [5]  and  others  for  colocated  schemes  that  solved  the  equations  in  an 
uncoupled  (pressure-correction)  manner.  However,  Chidambaram  [4]  appears  to  be  the 
first  to  derive  and  employ  an  appropriate  interpolation  scheme  for  the  compressible 
Navier-Stokes  equations  solved  as  a  coupled  system. 

An  alternative  promising  procedure  that  is  being  evaluated  at  present  employs  a  sixth 
order  compact  filter  [6]  that  appears  to  minimize  or  eliminate  grid-frequency  oscillations 
without  damping  the  solution  otherwise.  An  example  of  grid-frequency  oscillations  can 
be  seen  in  the  contours  of  mean  streamwise  velocity  in  a  rib-roughened  channel  shown  in 
Fig.  1.  Ribbed  channels  are  somewhat  characteristic  of  the  internal  cooling  passages  of 
gas  turbine  engine  blades.  The  LES  solution  employed  a  dynamic  subgrid-scale  model 


Figure  1.  Contours  of  streamwise  component  of  mean  velocity 
for  flow  in  a  rib-roughened  channel,  unfiltered 
solution. 

on  a  grid  composed  on  40  x  32  x  24  control  volumes.  The  simulation  was  at  a  Reynolds 
number  of  20,000  based  on  the  hydraulic  diameter.  The  ratio  of  the  rib  height  to  the 
channel  height  was  0.2  and  the  rib  streamwise  spacing  was  7.2  channel  heights.  Figure  2 
shows  the  velocity  contours  when  the  sixth  order  filter  was  employed.  Figure  3  compares 
the  present  LES  results  (with  filtering)  for  the  mean  streamwise  velocity  midway  between 
ribs  with  experimental  data  [7,8]  and  the  LES  simulations  of  Ciofalo  and  Collins  [9]. 

The  agreement  with  experimental  data  is  encouraging.  Notice  that  a  small  region  of 
recirculation  is  present  at  this  location,  3.6  rib-heights  downstream.  Figure  4  compares 
the  rms  streamwise  velocity  distribution  at  the  same  location.  Again,  the  agreement  with 
measurements  is  encouraging.  Although  not  shown  here,  the  effect  of  filtering  on  the 
mean  velocities  is  almost  indiscernible  on  plots.  Small  changes  can  be  seen  in  rms 
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Figure  4.  RMS  streamvvise  velocity  fluctuations 

midway  between  ribs  normalized  by  bulk 
velocity,  filtered  solution. 


quantities,  but  the  filtered  solution  tends  to  agree  best  with  the  experimental  data.  This  is 
work  in  progress,  and  full  results  have  not  yet  been  published. 

2.4  Incompressible  Flow  Results 

Although  the  main  thrust  of  the  research  involved  an  “all-speed”  compressible 
formulation,  some  aspects  of  the  approach  were  evaluated  first  through  computations 
with  an  incompressible  formulation  in  order  to  conserve  computer  resources.  Also,  it 
was  then  possible  to  too  take  advantage  of  the  greater  abundance  of  experimental  and 
computational  results  available  for  comparison  in  the  incompressible  regime. 
Furthermore,  before  attempting  LES  or  DNS  simulations  with  any  algorithm,  a  number  of 
traditional  two  and  three-dimensional  flows  were  computed  including  the  three- 
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dimensional  driven  cavity  and  the  growth  of  a  small  disturbance  in  a  two-dimensional 
channel.  Such  results  are  described  in  the  dissertations  by  Wang  [3]  and  Dailey  [10]. 
The  evolution  of  a  small  disturbance  was  considered  particularly  useful  in  screening 
candidate  discretization  schemes  for  temporal  and  spatial  accuracy. 

The  incompressible  planar  channel  flow  at  a  bulk  Reynolds  number  of  about  5500 
was  the  primary  flow  chosen  to  demonstrate  the  capability  of  the  schemes  developed  to 
perform  large  eddy  simulations  of  turbulent  flow.  The  merits  of  the  schemes  and 
subgrid-scale  models  were  gauged  by  comparisons  with  the  fine  grid  DNS  results  of  Kim 
et  al.  [11]  which  used  four-million  grid  points  and  the  experimental  work  of 
Niederschulte  et  al.  [12].  Eight  different  combinations  of  the  three  different  spatial 
algorithms  (i.e.,  2CD,  4CD  and  UPWIND),  combinations  of  the  regular/staggered  grid 
arrangements,  and  two  different  SGS  models  were  evaluated.  Effects  of  grid  refinement 
were  included  in  the  study.  Extensive  comparisons  between  the  calculated  results 
obtained  with  the  various  discretizations,  grids,  and  subgrid-scale  models  can  be  found  in 
[2,3]. 

A  second  configuration  simulated  with  the  dynamic  model  and  the  incompressible 
formulation  was  turbulent  flow  in  a  square  duct  at  a  Reynolds  number  of  about  9500 
based  on  the  bulk  velocity  and  duct  width.  To  the  principle  investigator’  knowledge,  this 
was  the  first  time  that  the  dynamic  model  had  been  used  for  this  configuration.  For  the 
square  duct  flow,  laboratory  experiments  have  indicated  that  if  the  flow  is  turbulent  the 
contours  of  mean  axial  velocity  bulge  toward  the  comer  due  to  a  secondary  flow  pattern 
in  the  cross-section.  This  phenomenon  is  not  observed  in  a  circular  duct  flow  or  in  the 
laminar  flow  of  a  non-circular  duct.  This  is  the  well-known  "secondary  flow  of  the 
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second  kind."  The  secondary  flow  is  caused  by  the  highly  anisotropic  and 
inhomogeneous  nature  of  the  turbulence  stresses  near  the  intersection  of  the  walls. 

Numerical  calculations  of  non-circular  duct  flows  using  the  Reynolds-averaged 
equations  were  first  performed  in  the  early  1970s.  However,  due  to  the  highly 
anisotropic  nature  of  the  turbulence  near  the  comer,  a  turbulence  closure  model  based  on 
the  flow  isotropy  assumption,  e.g,,  the  k  -  e  (two-equation)  model,  fails  to  resolve  the 
secondary  flow  (unless,  special  provisions  like  using  a  nonlinear  form  of  the  two- 
equation  model  [13]  was  employed).  A  second-order  Reynolds  stress  model  can 
recognize  the  flow  anisotropy.  Nevertheless,  the  empiricism  used  in  the  modeling  of  the 
various  correlation  terms  in  the  Reynolds  stress  transport  equation  has  restricted  the 
generality  of  the  model.  In  this  respect,  DNS/LES  may  provide  an  alternative  tool  to 
study  the  more  complicated  near-wall  turbulence  induced  by  the  secondary  flow. 

The  simulation  was  carried  out  on  a  65  x  65  x  65  grid  with  a  computation  domain  of 
2nD  X  D  X  D ,  where  D  is  the  duct  width.  Since  there  was  only  one  homogeneous 
direction,  sample  points  had  to  come  mainly  from  the  time  step  solution;  this  extended 
the  number  of  time  steps  needed  to  collect  enough  samples  to  make  the  statistics 
stationary. 

Figure  5  shows  the  cross-stream  secondary  velocity  vectors  in  the  lower  left 
quadrant.  Obviously,  the  secondary  flow  exhibits  a  high  degree  of  symmetry  with 
respect  to  the  comer  angle  bisector  (although  more  sample  points  are  probably  needed  to 
achieve  perfect  symmetry).  The  computed  secondary  flow  pattern  agrees  with 
experimental  observations  with  the  flow  being  toward  the  comer  along  the  comer  angle 
bisector  and  away  from  the  comer  along  the  walls.  The  maximum  amplitude  of  this 
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Figure  5.  Secondary  flow  velocity  vectors  in  the  lower-left 
quadrant,  square  duct,  Rcc  =  1 1,500,  dynamic 
model. 

secondary  velocity  was  1.5%  of  the  bulk  How.  Results  were  compared  with  data 
available  in  the  literature.  Additional  comparisons  on  this  case  can  be  found  in  [2,3]. 

Comparisons  were  also  made  with  incompressible  results  to  validate  the  several 
discretizations  and  models  employed  by  using  a  compressible  formulation  at  very  low 
Mach  numbers.  Flows  simulated  in  this  manner  included  the  decay  of  isotropic 
turbulence  as  well  as  the  channel  flow.  These  results  have  been  reported  in  [4,10,  14-16]. 


2.5  Results  for  Flows  with  Variable  Property  Heat  Transfer 
Two  heating/cooling  arrangements  have  been  studied  to  date  in  this  research.  The 
first  was  a  channel  heated  on  one  side  and  cooled  at  the  same  rate  on  the  opposite  side. 
Temperature  ratios  of  1.02  and  3.0  were  employed.  The  simulations  were  carried  out  for 
air  at  a  nominal  Mach  number  of  0.0 1  and  a  Reynolds  number  based  on  centerline 
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properties  and  channel  half  width  of  about  3200  and  a  Prandtl  number  of  0.71.  An 
implicit  compressible  upwind  scheme  that  was  fourth  order  accurate  in  the  viscous  terms 
and  third-order  accurate  in  the  advective  terms  was  used  on  -a  65  x  65  x  65  grid  that  spanned 
a  27i5x  25  X7t5  domain  where  5  is  the  channel  half  width.  The  grid  in  the  wall-normal 
direction  was  stretched  resulting  in  the  point  closest  to  the  wall  being  located  at  y""  =  0.92 . 
As  can  be  seen  in  Fig.  6,  the  van  Driest  transformation  did  not  collapse  the  simulation 
data  onto  the  incompressible  law  of  the  wall  for  the  case  with  the  highest  heat  transfer 
rate  (HIGH).  The  results  labeled  LOW  were  obtained  using  a  temperature  ratio  of 
1.02. 

Although  the  Mach  number  was  in  the  nominally  incompressible  range,  density  rms 
percentage  fluctuations  were  about  9%  near  the  cold  wall  (Fig.  7).  This  is  almost  as  large 
as  the  density  rms  fluctuations  observed  by  Coleman  et  al.  [17]  in  their  DNS  results  for 
channel  flow  at  a  Mach  number  of  3.0.  Although  this  was  initially  surprising,  it  becomes 
plausible  when  viewed  in  the  following  way.  It  is  reasonable  that  the  overall  temperature 
difference  in  the  flow  be  the  scaling  factor  to  determine  the  order  of  magnitude  of  the 
temperature  fluctuations.  Indeed,  when  the  distributions  for  the  two  cases  are 
nondimensionalized  by  half  the  difference  in  wall  temperatures  (curve  labeled  "original"), 
the  distributions  for  the  two  cases  do,  in  fact,  agree  in  order  of  magnitude  (but  not  detail) 
as  can  be  seen  in  Fig.  8.  That  being  the  case,  then  it  follows  that  on  a  local  percentage 
basis,  is  expected  to  achieve  a  larger  value  near  the  cold  side  than  near  the  hot  side. 
Note,  however,  that  on  an  absolute  basis,  the  peaks  near  the  hot  wall  would  be  about 
7%  larger  than  near  the  cold  wall.  Similarly,  the  velocity  rms  peaks  were  observed  to  be 
somewhat  larger  near  the  hot  wall  than  near  the  cold  wall.  Somewhat  unexpectedly,  the 


Kim  et  al. 

LOW  (cold) 

LOW  (hot) 
HIGH,  VD  (cold) 
HIGH  (cold) 
HIGH,  VD  (hot) 
HIGH  (hot) 


Figure  6.  Mean  velocity  profile  in  wall  coordinates  for 

channel  showing  effects  of  high  heat  transfer  rates. 
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Figure?.  Distribution  of  Pr^  across  channel;  scaled:  p^/(p); 
original:  P™,/p„. . 
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wall  shear  stress  was  observed  to  be  about  17%  greater  on  the  hot  wall  than  on  the  cold 
wall.  Other  relevant  results  such  as  the  shear  stress  and  heat  flux  budgets,  distributions  of 
~{TV) ,  skewness,  flatness,  and  correlation  coefficients  can  be  found  in  [3,18]. 

As  for  the  density,  this  study  as  well  as  others  at  higher  Mach  numbers  indicate  that 
the  mean  flow  is  nearly  isobaric.  If  it  is  assumed  to  be  so,  and  the  pressure  percentage 
fluctuations  are  assumed  to  be  small,  then  one  can  deduce  that  p^/{p}  =  T^I{T) .  By 
computing  the  density  rms  values  from  this  relationship,  one  obtains  agreement  within 
2%  of  the  results  shown  in  Fig.  7.  Coleman  et  al.  [17]  observed  similar  trends  for 
Pmi/(P)  and  T^I{T)  for  their  M  =  3  results,  i.e.,  an  indication  of  nearly  isobaric  behavior. 

Not  all  of  the  above  observations  were  expected,  so  a  follow-on  series  of  large  eddy 
simulations  were  conducted  in  order  to  shed  more  light  on  the  effects  of  property 
variations  on  turbulent  channel  flow.  For  this  purpose,  channel  flow  with  a  uniform  heat 
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flux  was  utilized  since  such  a  configuration  is  much  more  common  in  applications  than 
the  hot-cold  isothermal  wall  case  discussed  above.  The  uniform  heat  flux  flow  undergoes 
bulk  heating  or  cooling  so  that  periodic  boundary  conditions  in  the  streamwise  direction 
cannot  be  used  in  a  strict  sense.  In  the  present  simulations,  periodicity  in  the  streamwise 
density-velocity  product  was  assumed.  The  temperature  and  pressure  were  assumed  to  be 
step-periodic  with  the  temperature  step  being  computed  from  an  energy  balance  utilizing 
the  imposed  wall  heat  flux  and  a  mean  pressure  gradient  being  established  of  such  a 
magnitude  so  as  to  maintain  the  desired  mean  mass  flow  rate.  The  results  to  be  shown 
here  were  obtained  on  a  grid  that  used  48  x  64x  48  control  volumes.  The  value  of  at 
the  near-wall  control  volume  was  approximately  unity.  The  Smagorinsky  subgrid-scale 
model  was  used  (Cs  =  0.08)  with  van  Driest-type  damping.  Further  details  on  the 
numerical  formulation  can  be  found  in  [10,15]. 

Three  cases  have  been  computed,  all  for  an  inlet  Reynolds  number  of  approximately 
11,000  based  on  bulk  properties  and  the  hydraulic  diameter  and  a  Mach  number  of  0.001 
based  on  inlet  bulk  properties.  These  were  a  low  heating  =  1.025),  a  high  heating 

=  1.485 ),  and  a  high  cooling  (  T^/T^  =  0.564 )  case.  The  low  heating  was  selected 
in  order  to  compare  with  the  passive  scalar  DNS  results  of  Kasagi  et  al.  [19].  Most 
experimental  studies  of  heat  transfer  in  turbulent  internal  flows  have  been  carried  out  in 
tubes  rather  than  channels.  Kakac  [20]  points  out  that  noncircular  channel  friction  factor 
and  Nusselt  number  results  are  expected  to  be  within  10-15%  of  circular  tube  values  with 
the  parallel  plate  channel  results  being  7-11%  high.  Indeed,  that  trend  was  confirmed  by 
the  present  simulations.  The  Nusselt  numbers  based  on  bulk  properties  obtained  from  the 
simulations  were  34.52,  30.77  and  33.99  for  the  low  heating,  high  heating,  and  high 
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cooling  cases,  respectively.  The  corresponding  skin-friction  coefficients  were  0.007683, 
0.007822,  and  0.007787.  Comparisons  have  been  made  with  several  correlations  of 
experimental  data,  and  the  correlation  values  of  skin-friction  coefficient  and  Nusselt 
number  tended  to  be  within  ±8%,  20%,  and  13%  for  the  low  heating,  high  heating,  and 
high  cooling  cases,  respectively.  Generally,  the  Nusselt  numbers  based  on  bulk 
properties  from  the  simulations  did  not  show  as  much  sensitivity  to  wall  to  bulk 
temperature  ratios  as  indicated  by  the  correlations.  Most  of  the  experimental  data  were 
obtained  at  Reynolds  numbers  higher  than  used  in  the  simulations;  however,  several 
correlations  were  claimed  to  be  valid  for  the  Reynolds  number  of  the  present  simulations. 
Although  experimental  data  are  available  for  global  engineering  parameters  under  high 
heat  transfer  conditions,  relatively  little  information  exists  on  detailed  velocity  and 
temperature  profiles  and  the  statistical  features  of  such  flows. 

The  maximum  density  rms  percentage  fluctuations  were  again  observed  to  be  in  the 
9%  range  for  this  “incompressible,”  M  =  0.001  flow.  The  near-wall  rms  values  were 
greater  for  the  high  cooling  case  than  for  the  high-heating  case,  which  is  consistent  with 
the  observations  reported  above  for  the  hot-cold  isothermal  wall  channel  flow.  The 
distribution  of  the  temperature  rms  percentage  fluctuations  (not  shown)  was  nearly 
identical  to  that  of  the  density,  a  consequence  of  the  flow  being  essentially  locally 
isobaric.  One  of  the  most  relevant  findings  of  this  study  was  that  many  of  the  results 
from  the  simulations  for  both  heating  and  cooling  tend  to  collapse  toward  the  same 
distribution  when  plotted  in  semi-local  coordinates.  Such  a  collapse  had  been  noted  by 
others  for  high  speed  flow,  but  this  is  the  first  report  of  a  similar  trend  for  nominally 
“incompressible”  flow.  The  semi-local  coordinates  are  obtained  by  using  local  properties 
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(density  and  viscosity)  when  computing  the  friction  velocity  (wall  shear  stress  divided  by 
local  density)  and  y*.  An  asterisk  will  be  used  to  denote  the  semi-local  coordinates  rather 
than 

Figure  9  shows  the  mean  velocity  profile  in  conventional  wall  coordinates  for  the 
three  heat  flux  cases.  For  reference,  the  constant  property  DNS  data  of  Kim  et  al.  [11] 
and  the  experimental  data  of  Niederschulte  [12]  are  shown  on  the  figure.  A  clear 
dependence  on  the  wall  to  bulk  temperature  ratio  is  evident.  When  the  data  are  plotted  in 
semi-local  coordinates,  Fig.  10,  the  spread  in  the  curves  is  greatly  reduced.  A  similar 
effect  was  observed  for  the  mean  temperature  profile  and  again  the  data  were 
significantly  compressed  by  use  of  semi-local  coordinates.  The  remarkable  thing  is  that 
this  near  similarity  extends  to  distributions  of  fluctuating  quantities,  also.  Figure  11 
shows  the  distribution  of  rms  velocity  components  plotted  in  traditional  wall  coordinates, 
and  Fig.  12  shows  the  same  data  plotted  in  semi-local  coordinates.  Notice  how  the 
transformation  shifts  the  location  of  the  maximums  and  tends  to  collapse  the  three  rather 
different  curves  into  a  single  distribution.  Clearly,  the  collapse,  while  significant,  isn’t 
perfect.  Further  results  from  the  constant  heat  flux  simulations  are  described  in  [10,21]. 

2.6  Use  of  Unstructured  and  Zonal  Embedded  Grids 

Research  on  the  use  of  both  unstructured  and  zonal  embedded  grids  has  been 
completed  as  part  of  this  project.  Unstructured  grids  allow  the  flexibility  of  placing 
computational  cells  within  the  problem  domain  in  an  arbitrary  manner.  Typically,  cells 
composed  of  tetrahedra  or  hexahedra  are  used.  The  zonal  embedded  grids  considered 
here  employed  zones  of  hexahedral  cells.  Each  zone  was  treated  in  a  structured  manner 
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Figure  11.  Velocity  fluctuations  in  wall 
coordinates  with  uniform 
heat  flux  heating  and  cooling. 


Figure  12.  Velocity  fluctuations  in  semi¬ 
local  coordinates  with  uniform 
heat  flux  heating  and  cooling. 


and  interpolation  was  used  at  zone  interfaces.  Fine  grid  zones  can  be  placed  in  regions 
where  greater  spatial  resolution  is  required. 

The  use  of  such  grids  is  an  attractive  alternative  for  dealing  with  geometrically 
complex  flows.  As  computing  hardware  continues  to  improve,  DNS  and  LES  will  be 
applied  to  flows  of  increasing  complexity.  Thus,  it  was  considered  appropriate  to  begin 
to  understand  the  advantages  and  limitations  associated  with  the  use  of  unstructured  and 
zonal  embedded  grids  for  LES  and  DNS.  Also,  the  use  of  such  grids  provides  the  ability 
to  cluster  computational  cells  near  all  solid  boundaries  while  permitting  a  more  sparse 
distribution  in  the  interior  of  the  domain.  Then,  if  algorithm  accuracy  can  be  maintained 
on  the  unstructured  and  embedded  grids,  less  computational  effort  will  be  required  to 
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perform  a  simulation  of  given  accuracy  than  would  be  required  on  a  more  traditional 
structured  grid. 

The  unstructured  grid  scheme  has  been  applied  to  the  LES  of  isotropic  decaying 
turbulence  using  a  grid  composed  of  tetrahedrons  with  encouraging  results.  These  results 
were  reported  in  [14].  More  recent  unstructured  grid  results  will  be  reported  in  [16].  For 
flows  along  solid  boundaries,  the  use  of  zonal  embedded  hexahedral  control  volumes  has 
shown  much  promise  as  a  means  of  economically  achieving  simulations  at  high  Reynolds 
numbers.  The  scheme  allows  for  small  volumes  near  solid  boundaries,  but  coarsening  in 
all  three  directions  is  possible,  leading  to  good  accuracy  with  a  relatively  small  number  of 
cells.  It  is  believed  that  a  major  advance  in  the  simulation  of  high  Reynolds  number 
flows  is  in  progress.  The  only  other  known  LES  application  of  zonal  embedded  grids 
was  reported  by  Kravchenko  et  al.  [22].  They  solved  the  incompressible  equations  using 
a  spectral  method  in  two  directions  and  a  Galerkin  method  with  B-spline  basis  functions 
in  the  third  direction.  The  present  form  of  their  method  can  only  be  applied  to  a  very 
limited  class  of  problems  compared  to  the  present  finite  volume  method  for  solving  the 
fully  compressible  equations. 

Sample  LES  simulations  have  been  carried  out  for  channel  flow  using  a  zonal 
embedded  grid.  A  projection  of  such  a  grid  on  the  x-y  plane  is  shown  in  Fig.  13.  The 
grid  shown  was  used  for  a  flow  at  an  approximate  Reynolds  number  of  3,500.  There 
were  5  zones  across  the  channel,  symmetrically  positioned  about  the  channel  centerline 
containing  a  total  of  128,736  cells.  The  resolution  of  the  near  wall  region  is  64  x  9  x  64 . 
Domain  decomposition  can  be  easily  used  to  achieve  good  load  balance  on  parallel 
processing  machines  such  as  the  Origin  2000.  A  simple  structured  grid  with  the 
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Figure  13  Two-dimensional  projection  of  a 
zonal  embedded  grid. 

equivalent  near  wall  resolution  would  require  roughly  twice  as  many  grid  points  (64^)  as 
the  zonal  embedded  grid.  This  represents  considerable  savings  in  CPU  cost  (an  estimated 
factor  of  4-6). 

One  of  the  challenges  with  zonal  grids  is  to  provide  smooth  and  conservative 
transitions  across  zonal  boundaries.  Figure  14  shows  rms  fluctuation  results  for  the  Re  = 
3,500  case  mentioned  above.  Zonal  boundaries  are  marked  on  the  figure.  The  solutions 
appear  to  be  smooth  across  the  boundaries  confirming  the  merits  of  the  interpolation 
scheme  used  at  the  zonal  interfaces.  Simulation  results  obtained  with  both  the 
Smagorinsky  and  dynamic  subgrid-scale  models  are  shown  along  with  the  fine  grid  DNS 
results  of  Kim  et  al.  [11].  Further  results  using  the  zonal  embedded  grid  are  shown  in 
Figs.  15  and  16  for  channel  flow  at  the  much  higher  Reynolds  number  of  21,000  based  on 
the  mean  streamwise  velocity  and  channel  half-height.  Results  were  obtained  with  both 
the  Smagorinsky  and  the  dynamic  subgrid-scale  models.  Comparisons  are  made  with  the 
experimental  data  of  Wei  and  Willmarth  [23]  and  the  LES  results  of  Piomelli  [24]  and 
Najjar  and  Tafti  [25].  The  intermmediate  grid  results  were  obtained  using  a  near  wall 
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Figure  14.  Velocity  ims  fluctuations  (u,  v,  w 

components  from  top  to  bottom)  for  Re  = 
3,500,  zonal  embedded  grids. 


Figure  15.  Mean  stream  wise  velocity  in  wall 

coordinates,  Re  =  21,000,  zonal 
embedded  grids. 
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Figure  16.  Velocity  rms  fluctuations  for  Re  = 

21,000,  zonal  embedded  grids. 

resolution  of  70  x  70,  a  middle  zone  of  60  x  60,  and  a  core  region  of  50  x  50.  The  fine 
grid  employed  a  near  wall  resolution  of  80  x  80,  a  middle  zone  of  64  x  64,  and  a  core 
region  of  50  x  50.  Although  this  work  is  not  yet  complete,  it  appears  that  the  use  of 
zonal  embedded  grids  allows  greater  accuracy  to  be  achieved  for  a  fixed  number  of  cells 
than  is  possible  using  finite  difference  methods  with  traditional  grid  procedures. 
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